function condition_number()
% CONDITION_NUMBER 条件数计算和分析
% 
% 演示矩阵条件数的计算和其对线性方程组求解的影响

fprintf('=== 条件数分析 ===\n\n');

%% 1. 条件数的定义和计算
fprintf('1. 条件数的定义和计算\n');
fprintf('条件数 κ(A) = ||A|| * ||A^(-1)||\n\n');

% 良态矩阵示例
A_good = [4, 1; 1, 3];
fprintf('良态矩阵 A =\n');
disp(A_good);

cond_2 = cond(A_good, 2);      % 2-条件数
cond_1 = cond(A_good, 1);      % 1-条件数
cond_inf = cond(A_good, inf);  % ∞-条件数

fprintf('2-条件数: κ_2(A) = %.6f\n', cond_2);
fprintf('1-条件数: κ_1(A) = %.6f\n', cond_1);
fprintf('∞-条件数: κ_∞(A) = %.6f\n', cond_inf);

% 病态矩阵示例
fprintf('\n病态矩阵示例 - 接近奇异的矩阵:\n');
A_bad = [1, 1; 1, 1.0001];
fprintf('A =\n');
disp(A_bad);
cond_bad = cond(A_bad);
fprintf('条件数: κ(A) = %.2e\n', cond_bad);

%% 2. 条件数对线性方程组求解的影响
fprintf('\n2. 条件数对线性方程组求解的影响\n');

% 使用良态矩阵
b1 = [5; 4];
x1_exact = A_good \ b1;
fprintf('良态系统 A*x = b:\n');
fprintf('精确解: x = [%.6f; %.6f]\n', x1_exact);

% 给右端项添加小扰动
delta_b = [0.01; 0.01];
b1_perturbed = b1 + delta_b;
x1_perturbed = A_good \ b1_perturbed;

rel_error_b = norm(delta_b) / norm(b1);
rel_error_x = norm(x1_perturbed - x1_exact) / norm(x1_exact);

fprintf('右端扰动: δb = [%.2f; %.2f]\n', delta_b);
fprintf('扰动解: x = [%.6f; %.6f]\n', x1_perturbed);
fprintf('右端相对误差: %.6f\n', rel_error_b);
fprintf('解的相对误差: %.6f\n', rel_error_x);
fprintf('误差放大因子: %.6f\n', rel_error_x / rel_error_b);
fprintf('理论上界: κ(A) = %.6f\n', cond_2);

% 使用病态矩阵
fprintf('\n病态系统:\n');
b2 = [2; 2.0001];
x2_exact = A_bad \ b2;
fprintf('精确解: x = [%.6f; %.6f]\n', x2_exact);

% 给右端项添加相同的小扰动
b2_perturbed = b2 + delta_b;
x2_perturbed = A_bad \ b2_perturbed;

rel_error_b2 = norm(delta_b) / norm(b2);
rel_error_x2 = norm(x2_perturbed - x2_exact) / norm(x2_exact);

fprintf('扰动解: x = [%.6f; %.6f]\n', x2_perturbed);
fprintf('右端相对误差: %.6f\n', rel_error_b2);
fprintf('解的相对误差: %.6f\n', rel_error_x2);
fprintf('误差放大因子: %.2e\n', rel_error_x2 / rel_error_b2);
fprintf('理论上界: κ(A) = %.2e\n', cond_bad);

%% 3. Hilbert矩阵的条件数增长
fprintf('\n3. Hilbert矩阵的条件数增长\n');
fprintf('维数    条件数\n');
fprintf('----    --------\n');

dimensions = 3:10;
cond_numbers = zeros(size(dimensions));

for i = 1:length(dimensions)
    n = dimensions(i);
    H = hilb(n);
    cond_numbers(i) = cond(H);
    fprintf('%2d      %.2e\n', n, cond_numbers(i));
end

% 绘制条件数增长图
figure('Name', 'Hilbert矩阵条件数增长');
semilogy(dimensions, cond_numbers, 'b-o', 'LineWidth', 2, 'MarkerSize', 8);
xlabel('矩阵维数');
ylabel('条件数 (对数尺度)');
title('Hilbert矩阵的条件数增长');
grid on;

%% 4. 随机矩阵的条件数分布
fprintf('\n4. 随机矩阵的条件数分布\n');
n = 10;
num_trials = 1000;
cond_values = zeros(num_trials, 1);

fprintf('生成 %d 个 %d×%d 随机矩阵...\n', num_trials, n, n);

for i = 1:num_trials
    A_rand = randn(n, n);
    cond_values(i) = cond(A_rand);
end

mean_cond = mean(cond_values);
median_cond = median(cond_values);
max_cond = max(cond_values);
min_cond = min(cond_values);

fprintf('条件数统计:\n');
fprintf('平均值: %.2e\n', mean_cond);
fprintf('中位数: %.2e\n', median_cond);
fprintf('最大值: %.2e\n', max_cond);
fprintf('最小值: %.2e\n', min_cond);

% 绘制条件数分布直方图
figure('Name', '随机矩阵条件数分布');
histogram(log10(cond_values), 30, 'Normalization', 'probability');
xlabel('log_{10}(条件数)');
ylabel('概率密度');
title(sprintf('%d×%d 随机矩阵条件数分布', n, n));
grid on;

%% 5. 条件数估计
fprintf('\n5. 条件数估计\n');
fprintf('对于大型矩阵，精确计算条件数代价昂贵\n');
fprintf('MATLAB提供了条件数估计函数 condest\n\n');

% 创建一个较大的矩阵
n_large = 100;
A_large = randn(n_large) + 10 * eye(n_large);  % 对角占优矩阵

tic;
cond_exact = cond(A_large);
time_exact = toc;

tic;
cond_est = condest(A_large);
time_est = toc;

fprintf('矩阵大小: %d×%d\n', n_large, n_large);
fprintf('精确条件数: %.6f (耗时: %.4f秒)\n', cond_exact, time_exact);
fprintf('估计条件数: %.6f (耗时: %.4f秒)\n', cond_est, time_est);
fprintf('相对误差: %.2f%%\n', abs(cond_est - cond_exact) / cond_exact * 100);
fprintf('加速比: %.1f倍\n', time_exact / time_est);

end